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Abstract. - The dynamics of spins on a random graph with ferromagnetic three-spin inter- 
actions is used to model the compaction of granular matter under a series of taps. Taps are 
modelled as the random flipping of a small fraction of the spins followed by a quench at zero 
temperature. We find that the density approached during a logarithmically slow compaction 
- the random-close-packing density - corresponds to a dynamical phase transition. We discuss 
the role of cascades of successive spin-flips in this model and link them with density-noise power 
fluctuations observed in recent experiments. 



Analogies with glasses and spin glasses have long |3] been made to explain the complex 
properties of granular media, but it was not until the seminal experiments of the Chicago 
group Q| on granular compaction were carried out that serious theoretical attempts were made 
to quantify such analogies. These fall into roughly two classes: lattice-based models |3|.|4||| 
in a finite-dimensional space (which in general do not admit analytic solutions) or mean field 
models |||l0| (where each site interacts with a large number of other sites). In this Letter we 
argue that models on random graphs may be used to describe those aspects of the behaviour of 
granular matter which depend on the finite connectivity of the (disordered) grains, while still 
remaining analytically accessible. In particular this approach allows us to relate phenomena 
in the compaction of granular matter, such as the random-close-packing density, which is 
reached asymptotically after many taps, to the structure of the free-energy landscape of spin 
models on random graphs. 

The concept of geometric frustration is thought to be a crucial feature of the behaviour 
of granular matter. It refers to the existence of voids in the granulate, which cannot be filled 
in by neighbouring particles due to geometric constraints, either on the mobility of these 
particles, or on their compatibility in shape or size with the void. A class of models of this 
feature restrict the number of occupied neighbouring sites of a particle on a lattice, cither 
dynamically || or statically || . More generally, and expressed in the language of interacting 
spins (e.g. spin up for the occupied site of a lattice gas model), geometrical frustration may be 
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written in terms of ferromagnetic nearest-neighbour interactions and anti-ferromagnetic next- 
nearest neighbour interactions [Q. Models on random graphs may be viewed as the Bethe- 
approximation to finite-dimensional models of this type. The common feature of these models 
when translated into spin-language is the existence of multi-spin interactions on plaquettes 
on a random graph Jl2| . 

In the following we consider the simplest possible case of a multi-spin model defined on 
the plaquettes of a random graph, a 3-spin Hamiltonian on a random graph where N binary 
spins Si = ±1 interact in triplets 



where the variable Cijk = 1 with i < j < k denotes the presence of a plaquette connecting 
sites i,j, k and C^k = denotes its absence. Choosing djk = 1(0) randomly with probability 
2c/N 2 (1 — 2c/N 2 ) results in a random graph, where the number of plaquettes connected to 
a site is distributed with a Poisson distribution of average c. The connection with granular 
compaction is made in accordance with Edwards' thermodynamic hypothesis We inter- 
pret the local contribution to the energy in different configurations of the spins as the volume 
occupied by grains in different local orientations. In modelling terms, its random structure is 
an obvious advantage in the context of the disordered nature of granular media: Additionally, 
the locally fluctuating connectivity may be thought of as modelling the range of coordination 
numbers of the grains |Q . 

This Hamiltonian has recently been studied on a random graph in the context of satis- 
fiability problems in combinatorial optimization and on a 2D triangular lattice Jl5|. It 
has a trivial ground state where all spins point up and all plaquettes are in the configura- 
tion + + + giving a contribution of —1 to the energy. Yet, locally, plaquettes of the type 

h, — I — , H (satisfied plaquettes) also give the same contribution, although one may 

not be able to cover the entire graph with these 4 types of plaquettes in equal proportions. 
In this case, which occurs for c > c c ~ 2.75 |Q most ground states will be ferromagnetic - 
corresponding to a state with long-range order and a crystalline state of the granular medium; 

reaching the state + + + from h, — I — , H requires two spin- flips and a crossing of 

an energetic barrier. 

We argue that the 3-spins ferromagnetic Hamiltonian is a suitable modcl-Hamiltonian for 
two reasons: The mechanism by which energy barrier has to be crossed in going from one 
metastable state to another aims to model the situation where a grain in order to reach a 
void in a granular assembly has to push apart temporarily other grains thus increasing the 
volume. This has recently been argued to be an important ingredient in models of granular 
compaction fl6[| . 

The feature of the model responsible for the slow dynamics, however, is the degeneracy of 
the four configurations of plaquettes with SiSjSk = 1 resulting in competition between satis- 
fying plaquettes (dense arrangements of particles) locally (all states with even parity may be 
used, resulting in a large entropy) and global order (only the + + + state may be used). When 
grains are shaken, they rearrange locally, but locally dense configurations can be mutually in- 
compatible: voids appear between densely packed clusters due to mutually incompatible grain 
orientations (geometric frustration) - leading to competition between locally dense packings 
and global crystalline order. 

To mimic the action of tapping leading to compaction we choose the following dynamics: 
Each tap consists of two phases. First, in the dilation phase, a small fraction (taken to be 
10 -4 ) of the sites is chosen at random and have their spins flipped. This corresponds to the 
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Fig. 1 - Compaction curve at connectivity c = 3 for a system of 10 4 spins (one spin is flipped at 
random per tap). The data stem from a single run with random initial conditions and the fit (dashed 
line) follows (g) with parameters p^ — .971, po = .840, D = 2.76, and r = 1510. The long-dashed 
line (top) indicates the approximate density 0.954 at which the dynamical transition occurs, the long- 
dashed line (bottom) indicates the approximate density 0.835 at which the fast dynamics stops, the 
single-particle relaxation threshold. 



rapid acceleration of the sample during tapping. Then, the system undergoes several (in our 
case 3) Monte- Carlo sweeps at zero temperature: Each sweep consists of choosing TV sites in 
series and flipping their spins with probability 1 if the flip lowers the energy, with probability 
1/2 if the flip leaves the energy unchanged. This phase corresponds to the relaxation process 
after the tap, where the particles reach a new locally stable configuration. 

This dynamics is a simplified version, suitable for spin models, of the tapping dynamics 
used in cooperative Monte Carlo simulations of sphere shaking jL7|. It has been recently 
introduced independently in the context of spins models of granular compaction by Dean and 
Lefevre jl8| (however their 2-spin models lack the competition between local and global order). 
In the context of combinatorial optimization it corresponds to the Walksat-algorithm, known 
to be one of the most efficient algorithms for the solution of satisfiability problems |ll| . 

An example of a single run of the system defined by (|l|) is shown in Figure |l|. We can 
identify three regimes of this dynamics, first a very fast increase of the density up to a 
density po, then the slow compaction regime which takes the density up to poo, and finally an 
asymptotic regime. 

In the first regime, which only lasts a few taps, all sites orient their spins in parallel 
with the local field acting on that site. This corresponds to a fast dynamics whereby si ngle 
particles locally find the orientation maximizing the density leading to the density of po 
We therefore term this the single-particle relaxation threshold (SPRT). 

This fast regime may be illustrated by considering a single site i connected to 2k\ other 
sites and subject to the local field hi — ^/^^2j k CijkSjSk- For random initial conditions, the 
values of U = hiSi are binomially distributed with a probability of y 2 (l/2) ki if ki — U 

is even and zero if it is odd. If < zero-temperature dynamics will flip this spin, turn 1^ 



( x )During the fast regime one may also model the mobility of the grains by a dynamics of the bonds subject 
to kinetic constraints; the resulting graph typically differs from those constructed according to (hi), however 
the compaction dynamics does not change qualitatively. 
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to —li and turn (fej ± Z»)/2 satisfied (dissatisfied) plaquettes connected to it into dissatisfied 
(satisfied) ones. This will cause the lj of fc, ± /, neighbouring sites to decrease (increase) by 
2. This dynamics stops when all sites have I > giving po = 1/(37V) ^\ Zj. 

Neglecting loops in the random graph, the fast dynamics may be modelled by a simple 
population dynamics of TV elements, each with a Poisson distributed value of ki and a value of 
U distributed according to the initial binomial distribution. At each step a randomly chosen 
element with negative U has its li inverted, and ki ± U randomly chosen elements have their 
values of I decreased (increased) by 2 until Zj > Vi. Running the population dynamics 
for TV = 10 4 at c = 3 we obtain p = 0.835(1) which is shown as a dotted line in Figure 
|l|. Note that this density is found to be much higher than that of a typical 'blocked' state 
with Zj > Q Vi which is found to be 0.49 |2(|(see also the discussion in plj)); despite the 
fact that these states are exponentially dominant, the non-equilibrium, non-ergodic nature 
of the fast dynamics is responsible for taking the system to an (atypical) blocked state of 
higher density. At the end of the fast dynamics, there are no sites with more frustrated than 
unfrustrated plaquettes connected to them. The SPRT density thus appears as the density 
which is reached dynamically by putting each particle into its locally optimal configuration, 
as has also been found in lattice-based models M and simulations of sphere packings Jll| , 0] , 
which incorporate both fast and slow dynamics. 

The second phase of the dynamics consists of removing some of the remaining frustrated 
plaquettes and gives a logarithmically slow compaction Q leading from density p to 
Poo. The compaction curve may be fitted to the well-known logarithmic law 

p(t) = Poo - (Poo - Po)/(l + 1/D ln(l + t/r)) , (2) 

which may also be written in the simple form 1 + t(p)/r = cxp {D p ~ p " } , implying that 
the dynamics becomes slow (logarithmic) as soon as the density reaches po- With increasing 
density, free-energy barriers rise up causing the dynamics to slow down according to (|^) . The 
point where the height of these barriers scales with the system size marks a breaking of the 
ergodicity of the dynamics, a break-up of the phase-space into a large number (scaling expo- 
nentially with the system size) of disconnected clusters, and a saturation of the compaction 
curve. For small driving amplitudes, we thus identify the asymptotic density (random close 
packing) with a dynamical phase transition ||lo|,[l^,[2^1. 

To support this picture we give a simple approximation for the density Poo at which the 
dynamical transition occurs. Using the replica-trick InZ = lim n _,o d n Z n [E6| and standard 
manipulations we obtain for the average of the n-th power of the partition-function of the 
Hamiltonian (^) averaged the ensemble of random graphs 

« Z"))=n^ 1 ^( ( ?)exp|-/V^ C ( ( T)ln( C ( ; )) (3) 

+c/3-c/3 ^2 c(w)c(f)c(CT)exp{/3^w a rV a } 

where c(a) is an order parameter function defined on the domain of the 2™ vectors a a = ±1. 
A simple variational ansatz |l4[ |2?], ^8) implementing one step of replica-symmetry breaking 
(RSB) 

n/m ( p 

~f#\ - TT J 1 

J dh b G A (h h )[2cosh(ph b )} 
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Fig. 2 - The Fourier-transform of the preceding 1024 taps is used to produce a power spectrum as 
a function of time. We plot the power in the first octave (frequency 1/(1024 taps) -2/(1024 taps)), 
second octave (frequency 2/(1024 taps) - 4/ (1024 taps)), etc. up to the fourth octave (bottom to top). 
One finds that the fluctuations of the power in the different frequency bands are strongly correlated 
across many octaves, and correspond to sudden changes in the density (shown as dotted line). 

where G&(h) is a Gaussian with zero mean and variance A, gives the free energy subject to 
the variational ansatz f([3) — extr^m/iC/^ A, m) with 

a £ t o a \ /^(/3VAz)[2cosh(/3VAz)]™- 1 sinh(/3yA^) 

p Mp, A, ra) = - -= 

V ' /£>2[2cosh(/3VAz)]™ 

-^^ln(/ L>z[2cosh(/?VAz)] m ) - c/(3m) ln( J J j D Zl Dz 2 Dz 3 

8 cosh(/?VAzi) cosh(/3\/Az2) cosh(/3VAz 3 ) cosh(/3) + 8 smh(/3VA^i) smh((3VAz 2 ) smh(/3VAz 3 ) sinh(/3) 

where D(z) denotes the Gaussian measure with zero mean and variance one. The dynamical 
transition occurs at a temperature where d(Pf(p,A,m))/dm evaluated at m = 1 develops 
a minimum at finite A pq| . The corresponding density is marked with a horizontal line in 
Figure and agrees well with the asymptotic density reached by the tapping dynamics. 

The third regime of the compaction process is reached only asymptotically and consists of 
fluctuations around the random-close-packing density. 

Figure |l| shows marked fluctuations around the logarithmic compaction law, which have 
been the subject of detailed experimental investigations p9[] , where marked correlations be- 
tween the Fourier components of different frequencies of these fluctuations were found. In 
order to compare the results of our model with experiment we follow p9| and use Fourier 
transforms to plot the power of the timeseries p(t) in different frequency bands against time 
in Figure ^. As in the experimental data, the fluctuation of the power in a given band shows 
marked 'bursts' after periods of calm. As a result of these bursts the fluctuations in the the 
power spectrum are correlated over a wide range of frequencies: In figure |^ we have plotted, 
as in p9[ , the average of the correlation matrix Cy as a function of octave separation |i — j\, 

where Cy := Mij •\J~jj i . .^i)(m 3 ^IJ P rODes the non-Gaussian components of the correlation of 
the noise power in the ith and jth octaves. My is defined as the covariances of noise power 
fluctuations (SO^Oj) / ^((SO,) 2 ) ((SO,) 2 } for i ^ j and M u := ((S0 3 ) 2 ) / J2 ket (P k ) 2 , where 
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Fig. 3 - The rescaled covariances of the power fluctuations as are plotted as a function of the octave 
separation for both the ferromagnetic 3-spin model (squares) and the parking- lot model (diamonds). 
The definition of these quantities is provided in the text, a high value of the rescaled covariance 
indicates strong correlations of the power-fluctuations of two given frequency bands. 



the average is over 5000 time steps in the asymptotic regime, SOi are the power fluctuation 
around the average in the ith octave and Pk, k € i is the power in the kth frequency bin in 
octave i. Figure [| also shows the corresponding results for the parking-lot model pofl , where 
particles of a certain size are desorbed from and absorbed by a surface at random, subject to 
the constraint that none of the particles overlap. 

The experimental data analyzed in |2^] shows a slow decay of the average Cy with octave 
separation (particularly in the bottom and top of the sample) with decaying to about 0.6 
over 6 decades. This behaviour is reproduced by the 3-spin model, but not by the parking- 
lot model. We argue that the 'burst' found in the time-series of the density is resposible 
for the correlation of noise power over a wide range of frequencies. In spin models with 
finite connectivity such bursts arise quite naturally due to cascades of spin-flips. The crucial 
mechanism is that the flipping of a single spin alters the local fields acting on its neighbouring 
sites. The configuration of the spins on these sites may then no longer be locally stable, 
causing them to flip in turn. The first spin thus acts as a 'plug' releasing the neighbouring 
spins and setting off a cascade of successive spin-flips. A plug may also be composed of two 
or more sites, which need to have their spins flipped before neighbouring spins are released. 

Note that the cascade mechanism is entirely absent not only from generic fully connected 
models (where each spin interacts with all spins in the system, but with an interaction energy 
scaling as 1/yN) but also from the parking-lot model: Clearly the creation and filling of such 
a gap by a particle does not cause the appearance of further gaps, in the way spins flips may 
trigger a cascade. For this reason, the bursts observed in experiment and in our model are 
absent in the parking-lot model as noted in p9|| . In other, lattice-based models [|[^] of finitely 
connected granular media, however, cascades are also to be expected. 

To conclude, we have presented a finitely connected spin model of vibrated granular matter, 
which while reproducing the slow logarithmic relaxation associated with compaction, has also 
thrown new light on the fast relaxation mechanisms in granular dynamics. Within this model 
we have also identified the asymptotic density reached with a dynamical phase transition, and 
interpreted the density-noise power fluctuations observed in recent experiments in terms of 
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spin flip cascades. 
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